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Extracting reliable indications of chaos from a single experimental time series is a challenging task, 
in particular, for systems with many degrees of freedom. The techniques available for this purpose 
often require unachievably long time series. In this paper, we propose a new method to discriminate 
chaotic from multi-periodic integrable motion for a short time series provided it is very accurately 
measured. The method is based on analyzing higher order time derivatives of the time series. It 
exploits the fact that power spectra of chaotic time series exhibit exponential high-frequency tails, 
while, in the integrable systems, the power spectra are normally terminated at a finite frequency. 
We apply the above method to analysing signals generated by integrable and non-integrable systems 
of many interacting classical spins. 
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The problem of detecting deterministic chaos in an ex- 
perimental time series is of fundamental interest for the 
foundations of statistical physics [H It also arises in 
other contexts^!, e.g., in biomedical applications Q. In 
a typical difHcult situation, a system having many de- 
grees of freedom is suspected to be chaotic, but one can 
access the time evolution of only one degree of freedom. 
This means that the primary indicator of chaos, namely, 
the Lyapunov instability in the many-dimensional phase 
space cannot be investigated directly. A notable exam- 
ple in this regard was an attempt of Ref.Q to identify 
microscopic chaos in the measured trace of a Brownian 
particle. The approach of Ref. Q was to analyze the 
rate of information entropy (IE) production by this trace. 
The limiting value of this rate in the chaotic systems is 
known to be equal to the sum of the positive Lyapunov 
exponents. The results of Ref.Q were consistent with 
the possible presence of microscopic chaos, but, at the 
same time, were criticized as leaving open the possibil- 
ity that the same signatures might be produced by non- 
chaotic systems d, H]- Detection of microscopic chaos 
means that it should be discriminated from (i) a stochas- 
tic noise process characterized by the infinite limiting rate 
of IE-production and from (ii) a multi-periodic integrable 
motion characterized by the zero rate of IE-production. 
The difficulty here is that extracting the true limiting 
behavior of the rate of IE-production in many-particle 
systems typically requires unachievably long time series. 

The present paper focuses on issue (ii) above and pro- 
poses a new method to discriminate chaos from a multi- 
periodic non-chaotic motion in a short very accurately 
measured time series. We consider time series gener- 
ated by many-dimensional Hamiltonian systems, where 
both chaotic and non-chaotic motions are smooth in time, 
and where the non-chaotic motion is characterized by 
a sufficiently large number of frequencies that cannot 
be resolved by the Fourier transform of the time series. 
The method is based on analyzing the higher order time 



derivatives of the time series. It exploits the fact, which 
we demonstrate numerically, that the power spectra of 
time series generated by many-dimensional Hamiltonian 
chaotic systems have exponential high-frequency tails. 
To the best of our knowledge, the existence of these tails 
has never been proven rigorously, but otherwise reported 
in the studies of non-Hamiltonian or low-dimensional 
chaotic systems jil-fioj. In contrast, the power spectra of 
integrable systems normally have a high-frequency cut- 
off, and their high-frequency behavior is non-universal. 
Taking time derivatives of a time series progressively sup- 
presses the low frequency parts of the power spectra and 
amplifies the high-frequency part, thereby producing a 
number of clearly detectable signatures of chaos. Partic- 
ularly intriguing is the possibility to apply this method 
to the time series obtained by monitoring large quantum 
systems. Should a quantum time series produce the same 
signatures of chaos as expected for classical systems, such 
a result would shed a new light on the notion of quantum 
chaos. 

The above method is likely to be limited by the exper- 
imental noise (lo|. This limitation can be partially dealt 
with by using ultra-low noise sensing devices and by an 
appropriate high-frequency filtering of the signal. But, 
more importantly, it can be overcome by applying the 
method to the time series representing equilibrium fluc- 
tuations of extensive macroscopic variables, such as, e.g., 
the total magnetization of a spin system, in which case, 
the signal-to-noise ratio can be decreased by measuring 
these fluctuations in larger systems. 

We introduce our method by applying it to the time 
series of the total spin polarization for large clusters of 
interacting classical spins. These clusters are known to 
exhibit relaxati on prop erties [ll| that were linked to mi- 
croscopic chaos |12l4l4l |. but the interaction can also be 
chosen to lead to a non-chaotic multi-periodic dynam- 
ics. We mostly consider two 6x6x6 cubic clusters with 
periodic boundary conditions. One of these clusters is 
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integrable, while the other is not. We first demonstrate 
the chaotic character of the non-integrable cluster by cal- 
culating its largest Lyapunov exponent. After that, we 
generate rather long time series for the two clusters and 
demonstrate that one cannot discriminate between the 
chaotic and non-chaotic time series on the basis of the 
numerically accessible rate of IE-production. We then 
proceed with calculating the higher order time deriva- 
tives for the two time series and show that their seventh 
derivatives already look qualitatively different to discrim- 
inate the chaotic from the non-chaotic time series simply 
by visual inspection. We further introduce several quan- 
titative criteria characterizing this difference. Finally, we 
demonstrate the utility of our method by applying it to 
very short time series. 

The 6x6x6 spin clusters are characterized by the 
nearest-neighbor (NN) interaction Hamiltonian 
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where {Six, Siy, Siz) are the components of the ith classi- 
cal spin normalized by the condition Sf^ + Sfy + Sf^ = 1, 
and Jx, Jy and Jz are the coupling constants. For the in- 
tegrable cluster, we select the Ising Hamiltonian charac- 
terized hy Jx = Jy = and = 1, while, for the chaotic 
cluster, we select Jx = —0.65, Jy = —0.3, Jz = 0.7. In 
the Ising case, the motion is integrable, because the z- 
component of each spin is a constant of motion, while 
the X- and the j/-components simply precess in the local 
fields created by the frozen z-components of the neigh- 
bors. Therefore, the non-chaotic time series to be com- 
puted below is characterized by the superposition of 216 
different frequencies. The characteristic time scale of the 
dynamics of the both clusters is made equal since, in the 
both cases, J^ + Jy + Jz = ^ ■ The equations of motion 
for the Hamiltonian ([1]) were solved numerically using the 
discretization routine of Ref. 111. This routine conserves 



the energy of the system exactly. The typical discretiza- 
tion time step At was 0.01. The discretization errors 
were controlled by looking at the effects of the change of 
At on the Lyapunov exponents and on the power spectra 
of the time series. The initial orientations of spins were 
chosen randomly. Therefore, the energy of the system 
was close to zero. 

We have computed the largest Lyapunov exponent for 
both clusters using the method of Ref. • The method 
consists of choosing small initial distance do between two 
phase space trajectories, letting them diverge during time 
T and then resetting this distance back to do along the 
displacement direction just before the reset, and so on, 
repeating the above manipulation many times. If the 
system is chaotic, the spread of the two trajectories is 
eventually controlled by the largest Lyapunov exponent, 
which can be calculated as the limiting value of the ex- 
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FIG. 1: (Color online) Numerical values of the maximum Lya- 
punov exponent Xmax for the chaotic (a) and the integrable 
(b) systems as a function of the reset time r (dots). The 
lines represent the fits by a constant in (a) and by function 



in (b). The insets illustrate the growth of the dis- 



where j is the reset 



tance between a pair of trajectories during 5 reset intervals. 



index, k is the total number of resets, and dj is the dis- 
tance between the two trajectories just before the jth 
reset. For a fixed value of r, integrable systems can also 
produce small but finite value of A,nax because of the 
polynomial spread of the trajectories. However, for the 
polynomial spread, Amax 1/t (up to a logarithmic prcf- 
actor), while for the exponential spread the value of Amax 
should not depend on r. The dependences of Amax on r 
for the two clusters are shown in Fig. [TJ For the Ising 
cluster, Amax approaches zero approximately as 1/r as 
expected for an integrable system. The chaotic cluster 
exhibits r-independent value Amax = 0.63. The insets of 
Figs[TJa) and (b) show that distance growth between two 
resets is exponential in the chaotic case, and slower than 
exponential in the integrable case. 

We now compute for each of the two clusters the time 
series of length T = 1000 for the a;-component of the total 
spin polarization, Alx (t) and try to discriminate between 
chaotic and non-chaotic time series. The two time series 
look very similar as illustrated in Figs.[3l(a) and (b). 

In an attempt to detect the difference between the two 
time series, we investigated their rates of the (e,T) in- 
formation entropy production [l[. For this purpose, we 
"coarse-grained" the time and the magnetization axes 
in steps of r and e, respectively, to obtain a newly dis- 
cretized version of the time series (a stream of symbols). 
The (e,r) Shannon information entropy for patterns of 
length N is given by Hsh{e,T,N) = -Y.Pr log 
where i is the pattern index, and Pi is the pattern 
probability. The Shannon (e, r) entropy per unit time 
is defined as \\-misi^oahsh{^,T,N), where hsh{^,T, N) = 
i [Hsh{£,T,N + 1) - Hsh{£,T,N)\. For a chaotic sys- 
tem, hsh{^,T, N) is expected to approach the constant 
value equal to the sum of the positive Lyapunov ex- 
ponents in the limit r^-oo,iV^oo,e^'0 and 
T — > 0. However, for many-dimensional systems, the 
above limit is, typically, impossible to reach in practice, 
because, as N increases or e decreases, any finite time 
series quickly becomes too short to fairly represent the 
statistics of all possible patterns of length N. Once this 
happens, each pattern, that occurs, occurs only once. 
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FIG. 2: (Color online) Comparison between IE rates 
hsh{i,T, N) produced by the chaotic system (solid lines) and 
by the integrable system (dashed lines) for N — 2 and differ- 
ent values of t. 

and, as a result, hsh{i,T, N) 0. We, nevertheless, 
calculated /is/i(e, t, 2) in order to check if there is any 
indication of chaos before the effect of the finite length 
of the time series sets in. The results presented in Fig. [2] 
are nearly identical for both chaotic and non-chaotic time 
series. Similar results were also obtained for the Cohen- 
Procaccia entropy [l^. 

We now demonstrate, that one can easily identify chaos 
in the time series of Mx{t) by looking at its derivatives of 
the n-th order, which we denote as Mi"^(t). In Figs.lSJc) 
and (d), we exemplify this statement by presenting the 

(7) 

time evolution of the 7th time derivative Mx (t) for 
the seemingly indistinguishable time series appearing in 
Figs.[3l[a) and (b). (Th e plots for the lower-order deriva- 
tives are given in [16|.) For the chaotic time series, 
Mx'^\t) fluctuates noticeably faster than Mx{t) and has 
a rather random appearance, while, for the non-chaotic 
time series, Mx'^\t) has the appearance of a slowly mod- 
ulated periodic signal with period of the order of the 
characteristic time of Mx{t). 

The clear difference between Figs. [3{c) and (d) can 
be understood from the fact that the power spectrum 
of the original time series, Pico), and the power spec- 
trum of the nth time derivative P'^")(cj) are related as 
P(")(w) = w^"P(w). Both P{uj) and P(^)(w) for each 
time series are presented in Figs. [3^0) and (f). In the non- 
integrable case, P{uj) has exponential tail e"'*''"', where 
7 is a constant. As a result, for sufficiently large n, 

We propose to use this dependence as a quantitative sig- 
nature of chaos. Figure \Ele) includes the fit to of 
the form cj^^e"'''!"!. The important aspect of this depen- 
dence is not that it becomes exponential at sufficiently 
large frequences, but that, before it becomes exponential, 
it has the universally shaped broad maximum shifting 
with n to increasingly high frequencies. On the contrary, 
the power spectrum for the integrable cluster is termi- 
nated sharply at a certain maximum frequency uj^ax- As 
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FIG. 3: (Color online) Time series of Mx(t) for 6 x 6 x 6 
spin lattices. Figures in the left column represent the non- 
integrable system: (a) fragment of the time series Mxit); (c) 
corresponding 7th time derivatives Mx \t); and (e) the power 
spectra P{uj) [blue] and P^^'(tj)[red]. Figures (b,d,f) in the 
right column represent the same quantities for the Ising sys- 
tem. The black dashed line in (e) depicts fit to P'^'(a;) of 
the form a oj^^e"''" where a is a fitting parameter and 7 
is the high-frequency exponential decay constant for P(i^). 
The power spectra were obtained as described in [l^. Note 
that the tails of the power spectra at a;/27r > 4 in (e) and 
u)/2-K > 0.6 in (f) are affected by the spectral leakage from 
lower frequencies due to the finite length of the time series. 

a result, the shape of a;^"P(a;) for sufficiently large n be- 
comes narrowly peaked around uimax- Therefore, uimax 
becomes the carrier frequency for the modulations ob- 
served in Fig. [3)^d), while the inverse width of this peak 
characterizes the modulation time scale. 

We further propose two related criteria of chaos. The 
first of them characterizes the root-mean-squared (RMS) 
values of M^\t), denoted as Mrms- In Fig. SJa), we 
plot the quantity P„ = Mrms/Mrms^^ as a function of 
n. For the non-chaotic time series, P„ exhibits satura- 
tion, while for the chaotic time series it increases nearly 
linearly without apparent limit. The second criterion 
looks at the evolution of w„ defined as the square root of 
the variance for the positive-o; part of P^"''(a;) [wq corre- 
sponds to P{i^)]- In Fig.UJa), we plot Wn = Vn/vQ. Here 
the difference is that, in the chaotic case, Wn asymptot- 
ically increases with n, while in the non-chaotic case it 
decreases, asymptotically approaching zero. 

Finally, the fact that the higher-order derivatives look 
more random for the chaotic time series should also have 
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FIG. 4: (Color online) Statistical characteristics of the time 
series of Mx"\t) for the chaotic and the integrable systems: 
(a) ratios R„ of the RMS value of Mi"' (t) to the RMS value 
of Mi""^'(t); (b) relative width W„ of the power spectrum 
P*"'(aj) with respect to the width of P(a;). 



quantifiable consequences in terms of the rate of the 
(e, T)-entropy Droductionflij. 

In the power spectrum of the non-chaotic time series 
shown in Fig. ^e) one can distinguish some discrete fre- 
quencies, which is an indication of integrability. The 
chaos indicators proposed above are mainly intended for 
the situations when this discreteness is not discernible 
due to either too large number of discrete frequencies or 
too short length of the time series. In Fig. [51 we illustrate 
the effectiveness of our method by applying it to a very 
short time series of length T = 10 produced by two 4x4 
square lattice spin clusters with the same nearest neigh- 
bor coupling coefficients as the previously considered two 
clusters. Detecting chaos from such a short time series us- 
ing just the power spectrum of Mx{t) is difficult, because 
it is contaminated by the oscillations and the power-law 
decays associated with the short length of the time series. 
However, taking the 7th derivative reveals the expected 
qualitative difference between the chaotic and integrable 
systems and the good quality fit of form ([2]). 

One potentially promising application of our method 
would be to analyze the equilibrium fluctuations of nu- 
clear spin polarization in solids [l7j|. These fluctuations 
were measured by several NMR groups [3 although 
not yet with high enough accuracy. The exponential 
high-frequency tail in the Fourier transform of the NMR 
free induction decay (2^ suggests that the signatures of 
chaos may be directly observable in this case. 

In conclusion, we have demonstrated that examining 
the behavior of higher-order time derivatives is an effec- 
tive way to distinguish between chaotic and integrable 
systems on the basis of a finite time series, provided the 
time series can be measured with a sufficient accuracy. 
For the time series with zero average, the criterion of 
chaos can be based on Eq.([2|) and its observable conse- 
quences. Wc arc grateful to H. Kantz, P.Gaspard and D. 
Shcpclyansky for discussions. 
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FIG. 5: (Color online) Time series of Mx{t) for 4 x 4 spin 
lattices. Same notations as in Fig. [3] Frames (a,b,c,d) now 
represent the entire time series used to obtain frames (e,f). 
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SUPPLEMENTARY MATERIAL 



A. Comparison between derivatives of the time series produced by the chaotic and the integrable systems 
In Fig. 3 of the main article, we presented the fragments of two time series Mx{t) of the total spin polarization for 

(7) 

the chaotic and the integrable 6x6x6 spin clusters together with their respective 7th derivatives, Mx (t). In Fig. SI 
below, we show the plots for all the derivatives Mx^\t) of the same two time series up to the 9th order. We used the 
standard finite-difference numerical procedure for computing these derivatives. Namely, if the discretized time series 
in a given order was |..., (ti, M^"^^ , ^i^+i, M^"^_^-^^ , ^^^+2, M^"^_^2) j then the next-order derivative was obtained 



as — ^f^^Tif-^^J 7 l^ii+ii — ''t^2-t+i^^ j ' '"J ■ "^^^ numerical derivatives beginning with the ninth in the 

integrable case and tenth in the chaotic case show signs of numerical noise associated with the accumulated rounding 
errors. In order to exclude this concern, the presentation in the main article was limited to the 7th derivative. 




Chaotic 



Integrable 




10 15 
Time 



3x10- 



-3x10'' 



3x 10' 
2x 10 
1x10 
( 

-1x10 
-2x 10'' 
-3x10^ 




3x10 



-3x10'* 





10 15 20 25 
Time 



7 



Chaotic 



Integrablc 




10 15 20 25 
Time 

FIG. SI: Plots of Mx{t) and its nine derivatives for the chaotic system (left) and the integrable system (right). Note: the 
9th derivative for the integrable case exhibits extrinsic noise due to numerical rounding errors. 



B. Rate of the Cohen-Procaccia (e, r)-entropy production for the chaotic and the integrable systems 

It was pointed out in the main article that (i) the results similar to those presented in Fig. 2 can be obtained on 
the basis of the Cohen-Procaccia (CP) entropy [A. Cohen and I. Procaccia, Phys. Rev. A 31, 1872 (1985)]; and (ii) 
the apparent impression that the 7th derivative of the time series in Fig. 3 looks more random for the chaotic system 
than for the integrable system is likely to have signatures in terms of the rate the (e, T)-entropy production. In this 
section, we present evidence supporting both of the above statements. 

In Fig. 2, we show the CP (e, r)-entropy plots corresponding to the same data set as those analysed in Fig. 2 of 
the main article with the help of the Shannon entropy. The magnetization axis is not discrctized for the calculation 
of the CP (e, T)-entropy. Instead, the distance between two patterns of length N is given by the maximum of the 
differences between each pair of time-discretized points. To calculate the CP (e, r)-entropy, Hcp{£,t, N), a group 
of i? = 100 reference patterns of length N is selected randomly, and the probability of each of them is obtained by 
counting the number of all other patterns occurring in the time series which are within distance s from that reference 
pattern. The CP (e, T)-entropy is given by 

Hcp{e,T,N) = -^Y.^ogP,. (3) 

{R} 

The limiting rate of the CP (e, T)-entropy production is defined as limjv-i-oo hcp{e, r, N), where 

hcp{e, T,N) = - [His, T,N + 1)- Hie, r, N)]. (4) 

T 

In Fig. S2, we used t = 0.1 (10 discretization time steps) for one pattern element and varied N from 10 to 19 in 
different plots. It is evident from Fig. S2 that one cannot distinguish integrable from chaotic dynamics on the basis 
of numerically accessible behavior of hcpie,T, N) because of the finite length (T = 1000) of our otherwise very long 
time scries. 
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FIG. S2: Comparison between hcp{e,T, N) with r = 0.1 and various values of A*' for the chaotic time series (solid lines) and 
the integrable time series (dashed lines). Curves reaching higher maximum values of hcp correspond to the smaller values of 
N. 
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With regard to the seventh time derivatives of the both time scries, we do notice some quaUtative difference between 
the chaotic and the integrable cases, when we compare the CP (e, r)-entropy for each of these derivatives with the CP 
(e, r)-entropy for the original time series. In Fig. S3, we present the average of hcp{e,T, N) over different N for the 
original time series and for the 7th time derivatives using the set of values of N given in Fig. S2. The plotted values 
of e for hcp{e,T, N) associated with the 7th derivative were rescaled by dividing the true values by the following 
factors: 1.13 • 10''' in the chaotic case and 1.31 • 10'^ in the integrable case. These factors were chosen such to make the 
spread between the maximum and the minimum values of Mx {t) equal to that of the corresponding (t) thereby, 
in particular, compensating for the different dimensions of Mx{t) and Mx (t). Despite the apparent drawbacks of the 
above A^-averaging procedure, it does indicate, that the maximum accessible values of hcp{e,T, N) for fixed r and 
N are reduced for the integrable system, while staying about the same for the chaotic system. This suggests that 
the derivatives of the time scries for the integrable system produce less information when the order of the derivatives 
increases, while, for the chaotic time series, this is not the case. 
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FIG. S3: The A''-averaged values of hcp{t,T, N) calculated for both Mx{t) [blue] and Afi^'(t) [red] as described in the text. 
The values of e for Afi^' (t) are rescaled by the factors given in the text. 



C. Details on the calculation of the power spectra in Fig. 3 of the main article 

The power spectra in Fig. 3 of the main article were obtained by calculating the absolute value of the Discrete 
Fourier Transform (DFT) of the respective time series multiplied by a smooth-shaped window {wi} to mitigate the 
spectral leakage from low frequencies to high frequencies due to the finite length of the time series. That is, if the 

original time series is , ...|, then the modified time series is ^t;, WiA/^"''^ , ...|. We used the ten- 

percent Tukcy window [F. Harris, Proceedings of the IEEE 66, 51 (1978); J. Tukey, Spectral analysis of time series^ 
(Wiley, 1967)] defined as: 
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1 1 — cos 


27ri 
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.O.lA^t. 





for 0<i< OmNt 



1 for 0.05iVt < i < 0.957Vt, 



0.5 O - cos 



27r(i - Nt) 



for OMNt <i< Nt. 



where Nt is the index of the last discretized time point in {to, ...,ti, ...,tN^}. This window smoothly suppresses the 
time series to zero at to a-nd tjVj to reduce the finite-time- length effects. 
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We did not use the above procedure for the power series presented in Fig. 5 of the main article, because, in that 
case, the time series was very short (T = 10), and, as a result, the window function Wi would contaminate the relevant 
part of the power spectrum. 



